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Abstract 

Spectral properties and transition to instability in neutral delay differential equations are 
investigated in the limit of large delay. An approximation of the upper boundary of stability 
is found and compared to an analytically derived exact stability boundary. The approximate 
and exact stability borders agree quite well for the large time delay, and the inclusion of a 
time-delayed velocity feedback improves this agreement for small delays. Theoretical results 
are complemented by a numerically computed spectrum of the corresponding characteristic 
equations. 

1 Introduction 

An important area of research in physics and engineering is control theory, and a recent monograph 
by Scholl and Schuster jjj gives a good overview of the developments in this field. From a dynamical 
systems perspective, one could consider stabilisation of unstable fixed points or unstable periodic 
orbits [2, 3, 4J. In fact, even when these unstable periodic orbits are embedded in a chaotic attractor, 
they can still be stabilised by weak external forces, as it has been first proposed in a seminal paper 
by Ott, Grebogi, and Yorke [5]. Since then, several other methods of controlling unstable motion 
have been proposed. One of them is a time-delayed feedback control proposed by Pygaras in [6], 
which can be easily implemented in a wide range of experiments and is non-invasive, i.e., it vanishes 
as soon as unstable motion becomes stable [IlElElQiiiniQaQSlCg This method 

utilizes a difference between a signal at the current time and the same signal at some time ago. The 
scheme can be improved by introducing multiple time delays into the control loop |19j . Further 
considerations of multiple delay control, also referred to as extended time-delay autosynchronisation 
can be found in [U ED 1221 E3 EU • 

From theoretical point of view, introduction of a time delay into the system leads to an infinite- 
dimensional phase space and transcendental characteristic equations [251126] . This adds a significant 
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difficulty to the stability and bifurcation analyses of such systems. Some analytical results on time- 
delayed feedback control can be found, for instance, in [27J ESI GHl [3Ql GU G2]- In the case of linear 
time-delayed systems with non-delayed highest derivative, one can use the Lambert function to 
find the solutions of the corresponding characteristic equation [6'6\ 134] . However, for the neutral 
equations, i.e. equations which have time delays in the highest derivative term, this approach 
fails, and other approaches should be used. Furthermore, since neutral delay differential equations 
(NDDEs) often possess discontinuities in their solution, the numerical treatment and bifurcation 
analyses of such equations are much more involved than those of regular delay differential equations 
(DDEs). For instance, the existing packages for bifurcation analysis of DDEs, such as DDE- 
BIFTOOL [35] and PDDE-CONT [36] currently are unable to perform continuation for neutral 
systems. 

This paper is devoted to the analytical and numerical analysis of a time-delayed system of 
neutral type. Such models arise in a variety of contexts, such as biological and population dynamics 
models, see, for example, [37 1 138 1 [39], Balanov et al. [40] , for instance, derived a neutral DDE as a 
model for torsional waves on a driven drill-string. Another example studied by Blakely & Corron 
|41j is a model of a chaotic transmission line oscillator, in which an NDDE was used to correctly 
reproduce experimental observations of fast chaotic dynamics. The system to be studied in this 
paper was first introduced by Kyrychko et al. |42j in the context of hybrid testing, where it proved 
to be a good physical model for the description of the effects of actuator delays. It is noteworthy 
that in hybrid testing, quite often one encounters significant time delays due to actuator response 
time. Furthermore, actuator delay strongly depends on the stiffness of the system, hence it can 
vary considerably even in different runs of the same experiment. For this reason, the actual values 
attained by the actuator time delay can be quite high |441 1431 146] . 

It is important to note that NDDEs are different from DDEs in that they may possess a 
continuous as well as a point spectrum and their stability properties are far from being completely 
understood. Here, we investigate two different kinds of time-delayed feedback, both of which arise 
naturally in experimental settings. The first of these includes time delay in the feedback force, 
while the second introduces a velocity feedback. To understand the stability properties of the 
system, we will analyse asymptotic behaviour of the eigenvalue spectrum and identify regions of 
(in)stability in terms of system's parameters. The validity of these results will be compared to the 
numerical solution of the corresponding characteristic equation. For the particular system under 
investigation it is possible to find the stability spectrum analytically, and therefore it serves as a 
perfect test model for which it is possible to compare exact and approximate stability boundaries. 
It will be shown that although the approximation may deviate quite significantly from the exact 
boundary for small time delays, it gives a good agreement for larger time delays. Therefore, this 
approximation can be used for systems described by neutral DDEs with large time delays, where 
it is impossible to find the stability boundary analytically. 

2 Stability analysis 

2.1 Delayed force 

Consider the following NDDE [12]: 

z(t) + 2Qz{t) + z(t) + pz(t - r) = 0, (1) 

where dot means differentiation with respect to time t, and r is the time delay. In the context 
of hybrid testing experiments on a pendulum-mass-spring-damper system, £ stands for a rescaled 



2 



damping rate, andp is the mass ratio. Introducing v(t) = z(t) and u(t) = v(t)+pv(t), this equation 
can be rewritten as a system of differential equations with a shift: 



z{t) = u(t) — pv(t — r), 

u(t) = -2C[u(t)-pv(t-r)]-z(t), (2) 
v(t) = u(t) — pv(t — r). 

With the initial data (z(0), «(0)) = (zo,^o) € M x M and v(s) = 4>(s) G C[— r, 0], this system can 
be first solved on < t < r interval, then on r < t < 2r and so on, provided the following sewing 
condition is satisfied: 0(0) = uq — p<j){— r). This condition ensures that there are no discontinuities 
in the solutions at t = kr, k E Z + . For arbitrary initial conditions the sewing condition does not 
hold, and leads to jumps in the derivative of the solution [43J. 

The equation ([I]) has a single steady state z* = 0. The stability of this steady state is determined 
by the real part of the complex roots A S C of the corresponding characteristic equation 

A 2 + 2CA + 1 + pA V Ar = 0. (3) 

As it was already mentioned in the introduction, the existing bifurcation packages for the analysis 
of delay equations, such as DDE-BIFTOOL |35| and PDDE-CONT [36], currently do not provide 
capabilities of calculating eigenvalues for NDDEs. One of the reasons for this lies in the so-called 
behavioral discontinuity, a feature unique to NDDEs as compared to standard DDEs. This refers to 
the fact that even when all characteristic roots are stable for r = 0, for r being small and positive 
infinitely many of these roots may have unbounded real parts. In other words, a small variation of 
the time delay leads to an infinitely large root variation [UJ SB] • 

Several methods based on linear multi-step approach and pseudospectral differentiation have 
been recently put forward which provide an efficient tool for computing the characteristic spectrum 
of NDDEs [HI [501 EJJ • We have used this method to compute the spectrum of Eq. ^ , which 
is shown in Fig. [lj It can be observed that for small time delays (Fig. [if a)) the steady state is 
stable, as all the eigenvalues are in the left half-plane. As time delay increases, a pair of complex 
conjugate eigenvalues crosses the imaginary axis, as demonstrated in Fig. [TJb) ) and (c), leading 
to an instability. As time delay increases further still, the unstable eigenvalues return to the left 
half-plane, thus restoring the stability. 

In the case when the mass ratio p in Eq. exceeds unity, the steady state is unstable for any 
positive time delay r. It is worth noting that for \p\ < 1, the steady state may undergo stability 
changes/switches as the time delay is varied. To understand the dynamics of the system in the 
neighborhood of these stability changes, one can use the framework of pseudocontinuous spectrum 
used by Yanchuk et al. \31\ [52] for the analysis of scaling behavior of eigenvalues for large time 
delays, who followed an earlier work of Lepri et al. [53J on scaling of the spectra. Following this 
approach, one can express the asymptotic approximation of the eigenvalues for large r as 

x = ^ + i (n + -cp)+o(\), (4) 
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where 7, f2, and <f> are real-valued quantities, which are associated with the real and imaginary part 
of the eigenvalue A, respectively. Substituting this representation into the characteristic equation 
Q, gives to the leading order in C(l/r): 



1 - ft 2 + 2i(tt - P n 2 e-^e-^e- inT = 0. (5) 
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Figure 1: (Color online) Spectrum of the characteristic equation ^ for different time delays: (a) 
r = 2.5, (b) r = 3.32, (c) r = 5, and (d) r = 7. Parameter values are: £ = 0.1, p = 0.3. The solid 
lines show the asymptotic pseudocontinuous spectrum given by Eq. ([7]). 



±1, ±2, ±3,.... in equation ([5j), we can simplify this equation 

0. 



By choosing Q = = 2im/T, n 
to 

i - n 2 + 2i(n - P n 2 e -^e~^ = o. (6) 

From (|4j) it follows that Re(A) ~ r f(Q)/r and Im(A) ~ £1 up to the leading order, and therefore the 
eigenvalues A accumulate in the complex plane along curves (7(0), O), with the real axis scaling as 
rRe(A). Solving equation Q gives an expression for the real part 7 of the eigenvalue as a function 
of the Hopf frequency 
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A steady state can lose its stability via a Hopf bifurcation, at which point the tip of curve 7(0) 
will cross the imaginary axis. If this happens, there will be an interval of frequencies Oi < $7 < O2, 
for which 7(0) > and 7(^1) = 7(^2) = 0. This instability can be prevented, provided the interval 
of unstable frequencies Oi < f2 < O2 lies inside the interval [f2 n °, O n ° +1 ] for some no [31] ■ Here, 
fii 5 2 are two positive roots of the equation 7(0) = 0, which can be found from Eq. (FtT) as 



O 2 
"1,2 
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1 - 2C 2 ± v / (l-2C 2 ) 2 -l+P 2 
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For further analytical progress, we expand this expression for small values of £, which gives 
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Figure 2: (Color online) Comparison of the approximate upper bound of stability according to 
Eq. (10) (dotted line) with an exact stability boundary for £ = 0.1 in the (r,p) plane. The 



grayscale (color code) encodes the value of the largest real part of the complex eigenvalues A. 
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Since the actual values of the frequencies are 0( n ) = 27rn/r for any integer n, the distance between 
any two successive frequencies is 2tt/t, and hence the necessary condition for stability Af2 < 2tt/t 
can be written as 

1 1 CV 1 1 \ 

Tt^p ~ VT+p ~ J \VT+p + V^p) < 27r/r ' (10) 

For large enough time delay r, p asymptotically approaches a lower bound of stability which 
corresponds to Ail = 0. It can be obtained from Eq. (Tsl) by using 







(^1-^2)^1 + ^2) 



n 2 2 
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which yields 
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2Ca/1 - C 2 ~ 2C- 



(11) 



Figure 2 shows the plot of the approximate stability boundary (10) as a function of time delay 



r for a given small value of the damping £. The grayscale (color code) in this figure indicates the 
value of the largest real part of the eigenvalues in the spectrum of the characteristic equation Q 
for each value of p and r. As it follows from Fig. [2j the analytically derived formula (10) for the 
maxima on the stability boundary deviates from the exact stability peaks (which correspond to 
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Figure 3: (Color online) Spectrum of the characteristic equation (13) or different time delays: (a 



t = 2, (b) r = 2.725, (c) r = 3.5, and (d) r = 5. Parameter values are: & = 0.25, (2 = 0.24, 



p = 0.3. The solid lines show the asymptotic pseudocontinuous spectrum given by Eq. (15) 



codimension-two Hopf bifurcation) for small delays, but it provides a good approximation for large 
time delay r. Appendix A contains an exact analytic expression for the stability boundary in terms 
of system's parameters. 



2.2 Delayed viscous damping 

In order to analyse the influence of velocity feedback on the stability of NDDE, we modify Eq. ([!]) 

as 

z(t) + 2&z(i) + z{t) + pz(t - r) + 2( 2 z(t - r) = 0. (12) 

This equation was introduced in Ref. [42] . where it was shown that depending on the difference 
between two damping parameters £i and (2, the stability domain may shrink and even split into 
separate stability regions in the parameter plane (the so-called death islands). The characteristic 
equation now modifies to 

A 2 + 2CiA + 1 + p\ 2 e- XT + 2C 2 Ae" Ar = 0. (13) 

Figure [3] shows the numerical approximation of the roots of this equation in the neighborhood of 
the origin. From this figure it follows that similar to the situation without velocity feedback, the 
system undergoes successive stability switches as the time delay is varied. 
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Figure 4: (Color online) Comparison of the approximate upper bound of stability according to 



Eq. ( 17) (dotted line) with an exact stability boundary for (j = 0.25 and (,2 = 0.24 in the (r,p) plane. 
The grayscale (color code) encodes the value of the largest real part of the complex eigenvalues A. 



Assuming in Eq. ( 13 ) the same asymptotic behavior Q of the eigenvalues for large time delay 
(i.e., the real part of the eigenvalue scales as 1/r), gives to the leading order 



0. 



(14) 



with the constraint VL = = 2im/T, n = ±1, ±2, ±3,.... One can solve this equation for the 
real part of the eigenvalue 7 at the Hopf bifurcation as a function of frequency O, as 



7 (fi) 



1 (l-ny + ^fo 2 

2 11 p^ + 4C 2 2 ^ 2 ' 



(15) 



Transition to instability occurs when 7(0) = 0, which gives the expression for instability fre- 
quencies 



O 2 



1 + 2(1 - 2(f ± 



l + 2C 2 -2C 2 ) 2 -l+p 



1 -p 2 

In a manner similar to the analysis of the delayed force feedback, one can make further analytical 
progress by assuming that both damping coefficients are small: \(i\ <C 1, IC2I *C 1. The necessary 
stability condition Af2 = f^i — f2 2 < 2-k/t gives the following asymptotic approximation for the 
maxima of the stability boundary: 
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< 2vr/r. 



(17) 
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The expression (17) can be further simplified for large time delay in a manner similar to ( 11 ), which 
2vW 



C|)(i-C 1 2 + C|)-2 v ^f 



t 2 



gives p 

It is noteworthy that the inequality (17) provides a good approximation for the stability bound- 
ary even when actual values of damping coefficients (j and £2 are large, as long as the difference 
(£ 2 — £f ) is small by the absolute value. Figure [4] shows an excellent agreement between the 
asymptotic approximation (17) and the exact stability boundary, especially for sufficiently large 
time delay. In Appendix A it is shown how the exact stability boundary varies depending on 
parameters, and in particular on the relation between the two damping coefficients. 



3 Conclusions 



Time delays are an intrinsic feature of many physical, biological and engineering systems, and 
in recent years the analysis of such systems has led to many interesting and important findings. 
There are systems where the time delay is present intrinsically due to processing times, mechanical 
inertia etc., and there are those where the time delay is introduced externally in order to stabilise 
unstable periodic orbits and steady states. Therefore, a better understanding of delay differential 
equations will provide a clear picture of the system's stability and controllability. In this paper we 
have concentrated on the analysis of two neutral delay differential equations. We have shown that 
depending on the time delay r, the systems exhibit stability switches, where stability is lost /regained 
depending on the time delay. In the case of delayed velocity feedback, the interplay between the 
time delay and the two damping coefficients gives different stability regimes in the parameter plane, 
and for some parameter values the stability area collapses into separate islands. We have derived 
an asymptotic approximation of the stability peaks for large time delays, based upon universal 
scaling arguments, and have compared this approximation with the exact stability boundary. The 
results agree quite well even when the time delay is not too large, and give excellent agreement for 
large delays. The results presented in this paper include numerical simulations of the characteristic 
spectrum and constitute the first attempt to approximate stability peaks for neutral DDEs. As it 
has already been mentioned in the introduction, neutral DDEs arise naturally in a wide range of 
physical problems, which makes the approach developed in this paper a useful tool for the stability 
analysis of such systems. 
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A Exact stability boundary 

To find an exact analytical expression for the stability boundary, one has to substitute A = icu into 
the characteristic equation ([3]). After separating real and imaginary parts, this gives [53] 

1 — uj 2 — puj 2 cos(wr) = 0, , , 

2(u + puJ 2 sin(ojT) = 0. ^ ' 

Squaring and adding these equations gives 

(l-p 2 )u; 4 + (4C 2 -2)w 2 + l = 0, (19) 
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Figure 5: Exact stability boundary of the characteristic equation Q in (r,p) parameter plane for 
£ = 0.1. The steady state is stable in the shaded area. 



This equation can be solved as 



W?2 



l-p 2 



l-2( 2 ± VU-2C 2 ) -1+P 2 



(20) 



In fact, Eq. (19) provides an expression for stability boundary value of p as parametrized by 



the Hopf frequency ui: 



= -^y/u* + 2u>* (2C 2 -1) + 1. 



The corresponding value of the time delay at the stability boundary is derived from Eq. (|18|) 



1 



2 C^ 

Arctan — r ± 7rA; 



1 



(21) 



(22) 



where k = 0, 1, 2, ... and Arctan denotes the principal value of arctan. Figure [5] illustrates the 
dependence of critical mass ratio p on the time delay r which ensures the stability of the steady 
state. It is noteworthy that if \p\ > 1, the steady state is unstable for any positive time delay r; 
on the other hand, if \p\ < 1 and £ > l/y/2, then the steady state is asymptotically stable for any 
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Figure 6: Exact stability boundary of the characteristic equation (13) in (r,p) parameter plane. 
The steady state is stable in the shaded area. Parameter values are: (a) £i = 0.25 and £2 = 0.24, 
(b) Ci = 0.23 and ( 2 = 0.25. 



positive time delay r. For ( < l/y/2, there is a lower bound on the value of p m i n = 2£y / l — ( 2 , so 
that an asymptotic stability is guaranteed for all r > provided p < p m i n [42J. 



In the case of time-delayed viscous damping, the characteristic equation ( 13 ) at the points of 
stability changes can be written as 



1 — lo 2 — pco 2 cos lot + 2&lo sin lot = 0, 
2(i + pcu sin lot + 2^2 cos lot = 0. 



(23) 



Squaring and adding these two equations gives the following parametrization of p by the Hopf 
frequency: 



"1,2 



1 



l-p 2 



1 - 2d 2 + 2C 2 2 ± V (1 " 2 Ci 2 + 2C 2 2 ) 2 " 1 +P 2 



(24) 



Similar to the previous case, one can derive parametric expressions for p(uS) and t(lo) from Eq. (23): 



p(co) 



LO' 



l) 2 + 4( Cl 2 -Cl 



The corresponding value of the time delay at the stability boundary can be found as 



t(lo) 



1 



LO 



2im — arccos 



P (l - LO 2 ) - 4ClC 2 

P 2 lo 2 + 4C| 



n = 1,2,3. 



(25) 



(26) 



Figure [6] demonstrates how stability boundary is affected by the relation between £1 and £2- In 
particular, we note that when (j = £2 the stability boundary touches the r-axis (p = 0), and for 
C2 > Ci) the stability area consists of non-overlaping death islands, inside which the oscillations are 
damped, and the steady state is stable. 
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